*******************************************************************
*
* 		Using community health workers to deliver a 
*	scalable integrated parenting program in rural China: 
*			a cluster-randomized controlled trial 
*
* Goal: Analysis for health paper to submit to SS&M
* Input: Luo et al. (2019) data.dta
* Author: Dorien Emmers
* Created on 04/07/2017
* Last modified on 09/10/2019
*
********************************************************************

clear all
clear mata
clear matrix
set maxvar 32000
set segmentsize 1g
set more off
macro drop all

* Set directories
global datadir 		""
global workingdir 	""
global outputdir 	""


* Import data
use 		"$datadir\Luo et al. (2019) data.dta", clear 

/*----------------------------------------------------------------------
Table 1: Baseline balance table
-----------------------------------------------------------------------*/
{		
orth_out 	cog_score2015 rc_score2015 ec_score2015 fm_score2015 gm_score2015 soemo_score2015 ill_diarrhea_bl ill_tract_bl ill_times_bl hb_bl ///
			ch_anemic_bl height_bl weight_bl waz_bl underweight_bl haz_bl whz_bl ///
			using  "$outputdir/tab1_balance", by(treatment) covariates(provid*) vce(cluster villid) ///
			pcompare count se numlabel armlabel("Control" "Treatment") title(Table 1: Baseline development, health, nutrition, and growth outcomes of intention-to-treat population)  replace
}
/*----------------------------------------------------------------------
Figure 2: ITT effects on Summary Indices
-----------------------------------------------------------------------*/
{
		eststo clear 

		eststo: reg knowledge_index		treatment 	knowledge_index_bl   		i.provid		, vce(cluster villid) 
		su knowledge_index if treatment==0
		estadd scalar cont_mean = `r(mean)'
		
		eststo: reg investment_index   treatment 	investment_index_bl    		i.provid		, vce(cluster villid) 
		su investment_index if treatment==0
		estadd scalar cont_mean = `r(mean)'
		
		eststo: reg food_index 			treatment	food_index_bl 				i.provid 		  , vce(cluster villid) 	
		su food_index if treatment==0
		estadd scalar cont_mean = `r(mean)'
		
		eststo: reg iron_rich_index		treatment	iron_rich_index_bl			i.provid 		 , vce(cluster villid) 	
		su iron_rich_index if treatment==0
		estadd scalar cont_mean = `r(mean)'
		
		eststo: reg meal_freq_index		treatment	meal_freq_index_bl			i.provid 		 , vce(cluster villid) 	
		su meal_freq_index if treatment==0
		estadd scalar cont_mean = `r(mean)'
		
		eststo: reg supplement_index	treatment	supplement_index_bl			i.provid 		, vce(cluster villid) 	
		su supplement_index if treatment==0
		estadd scalar cont_mean = `r(mean)'
		
		eststo: reg health_index		treatment	health_index_bl 			i.provid 		, vce(cluster villid) 	
		su health_index if treatment==0
		estadd scalar cont_mean = `r(mean)'
		
		eststo: reg skill_dev_index	treatment		skill_dev_index_bl			i.provid 	  i.tester_id_2016 , vce(cluster villid) 			
		su skill_dev_index if treatment==0
		estadd scalar cont_mean = `r(mean)'
		
		eststo: reg phys_growth_index	treatment	phys_growth_index_bl				i.provid		, vce(cluster villid)  
		su phys_growth_index if treatment==0
		estadd scalar cont_mean = `r(mean)'
			
						estout using  "$outputdir/tab_ITT_overview.xls", replace cells(b(fmt(2))  se(fmt(2))  ci(fmt(2)) ) starlevels(* 0.1 ** 0.05 *** 0.01)  mgroups("Knowledge index" "Investment index" "Food index" "Iron rich food index" "Meal frequency index" "Supplement index" "Health index" "Skill development index" "Physical growth index", pattern(1 1 1 1 1 1 1 1 1)) collabels(none) numbers("[" "]") mlabels(none) ///
						keep(treatment)	stats(pval1 r2 N cont_mean controls, labels("P-value" "R-squared" "No. of observations" "Control mean" "Controls"))											///
						title(ITT Effects on Summary Indices) 								///
						note("Expressed in z-scores. Z-scores were obtained by non-parametrically standardizing within age groups. Effects on primary outcomes are estimated intervention effects controlling for baseline scores. Standard errors are adjusted for clustering at the village level and stratification at the township level.")
		
		rwolf 			knowledge_index investment_index food_index iron_rich_index meal_freq_index supplement_index health_index skill_dev_index phys_growth_index, indepvar(treatment) bl(_bl) strata(provid) controls(i.provid) vce(cluster villid) reps(100) seed(563)
}									
						
																							
/*----------------------------------------------------
Table 3: ITT effects in Child Development
-----------------------------------------------------*/																			
{

* Panel A: Skill development
								
			eststo  clear 
			
			eststo: reg skill_dev_index				treatment	skill_dev_index_bl 		i.provid i.tester_id_2016 	, vce(cluster villid) 
			su skill_dev_index if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg cog_score		 			treatment	cog_score_bl					i.provid i.tester_id_2016 	, vce(cluster villid) 
			su cog_score if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg rc_score		  			treatment	rc_score_bl		 				i.provid i.tester_id_2016 	, vce(cluster villid)
			su rc_score if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg ec_score		  			treatment	ec_score_bl		 				i.provid i.tester_id_2016 	, vce(cluster villid)
			su ec_score if treatment==0
			estadd scalar cont_mean = `r(mean)'

			eststo: reg gm_score		  			treatment 	gm_score_bl	    				i.provid i.tester_id_2016 	, vce(cluster villid) 
			su gm_score if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg fm_score		   			treatment 	fm_score_bl		     			i.provid i.tester_id_2016 	, vce(cluster villid) 
			su fm_score if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg soemo_score			   		treatment	soemo_score_bl					i.provid  i.tester_id_2016	, vce(cluster villid)  
			su soemo_score if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			
							estout using  "$outputdir/tab_ITT_development.xls", replace cells(b(fmt(2))  se(fmt(2))  ci(fmt(2)) ) starlevels(* 0.1 ** 0.05 *** 0.01)  mgroups("Skill development index" "Cognition" "Receptive language" "Expressive language" "Gross motor" "Fine motor" "Social-emotional", pattern(1 1 1 1 1 1 1 1 1 1 1)) collabels(none) numbers("[" "]") mlabels(none) ///
							keep(treatment)	stats(pval1 r2 N cont_mean, labels("P-value" "R-squared" "No. of observations" "Control mean"))											///
							title(ITT Effects on Skill Development) 								///
							note("Expressed in z-scores. Z-scores were obtained by non-parametrically standardizing within age groups. Effects on primary outcomes are estimated intervention effects controlling for baseline scores. Standard errors are adjusted for clustering at the village level and stratification at the township level.")
								
					
			rwolf 			cog_score rc_score ec_score gm_score fm_score soemo_score, indepvar(treatment) strata(provid) controls(i.provid i.tester_id_2016) bl(_bl) vce(cluster villid) seed(563)	
			
	
* Panel B: Health

			eststo clear 
			
			eststo: reg health_index	 		treatment	health_index_bl 				i.provid 		, vce(cluster villid) 	
			su health_index if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg ill_diarrhea			treatment	ill_diarrhea_bl 				i.provid 		, vce(cluster villid)
			su ill_diarrhea if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg ill_tract	 			treatment 	ill_tract_bl					i.provid 		, vce(cluster villid) 
			su ill_tract if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg ill_times	 			treatment 	ill_times_bl	   				i.provid 		, vce(cluster villid) 
			su ill_times if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
							estout using  "$outputdir/tab_ITT_health.xls", replace cells(b(fmt(2))  se(fmt(2))  ci(fmt(2)) ) starlevels(* 0.1 ** 0.05 *** 0.01)  mgroups("Health index" "Diarrhea" "Respiratory Tract Infection" "No. of times ill" , pattern(1 1 1 1 1 1 1 1 1 1 1 1 1)) collabels(none) numbers("[" "]") mlabels(none) ///
							keep(treatment)	stats(pval1 r2 N cont_mean, labels("P-value" "R-squared" "No. of observations" "Control mean"))											///
							title(ITT Effects on Health) 								///
							note("Expressed in z-scores. Z-scores were obtained by non-parametrically standardizing within age groups. Effects on primary outcomes are estimated intervention effects controlling for baseline scores. Standard errors are adjusted for clustering at the village level and stratification at the township level.")
							
			
			rwolf 			ill_diarrhea ill_tract ill_times, indepvar(treatment) strata(provid) controls(i.provid ) bl(_bl) vce(cluster villid) seed(563) reps(100)
			
	
* Panel C: Nutrition and growth																	

			eststo  clear 
			
			eststo: reg hb		 				treatment	hb_bl							i.provid 	, vce(cluster villid) 	
			su hb if treatment==0
			estadd scalar cont_mean = `r(mean)'

			eststo: reg ch_anemic		 		treatment	ch_anemic_bl		 			i.provid 	, vce(cluster villid) 	
			su ch_anemic if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg phys_growth_index		treatment	phys_growth_index_bl			i.provid	, vce(cluster villid)  
			su phys_growth_index if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg haz		   				treatment	haz_bl 							i.provid	, vce(cluster villid)  
			su haz if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg healthy_haz		   		treatment	healthy_haz_bl 					i.provid	, vce(cluster villid)  
			su healthy_haz if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg healthy_whz				treatment	healthy_whz_bl 					i.provid	, vce(cluster villid)  
			su healthy_whz if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			eststo: reg healthy_waz				treatment	healthy_waz_bl 					i.provid	, vce(cluster villid)  
			su healthy_waz if treatment==0
			estadd scalar cont_mean = `r(mean)'
			
			
							estout using  "$outputdir/tab_ITT_nutrition_growth.xls", replace cells(b(fmt(2))  se(fmt(2))  ci(fmt(2)) ) starlevels(* 0.1 ** 0.05 *** 0.01)  mgroups("Hb" "Anemia" "Physical growth index" "HAZ" "Healthy HAZ" "Healthy WHZ" "Healthy WAZ", pattern(1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)) collabels(none) numbers("[" "]") mlabels(none) ///
							keep(treatment)	stats(pval1 r2 N cont_mean, labels("P-value" "R-squared" "No. of observations" "Control mean"))											///
							title(ITT Effects on Health) 								///
							note("Expressed in z-scores. Z-scores were obtained by non-parametrically standardizing within age groups. Effects on primary outcomes are estimated intervention effects controlling for baseline scores. Standard errors are adjusted for clustering at the village level and stratification at the township level.")
			
			rwolf 			haz healthy_haz healthy_whz healthy_waz hb ch_anemic, indepvar(treatment) controls(i.provid ) bl(_bl) vce(cluster villid) seed(563) reps(100)
						
}



/*----------------------------------------------------------------------
Table A.1: Baseline balance table of characteristics of ITT population
-----------------------------------------------------------------------*/
{		
orth_out 	agemonth sex2015 ch_premature minority ch_firstborn mom_age mom_at_home mom_edu hh_dibao hh_income ///
			using  "$outputdir/taba1_balance", by(treatment) covariates(provid*) vce(cluster villid) ///
			pcompare count se numlabel armlabel("Control" "Treatment") title(Table A.1: Baseline characteristics of intention-to-treat population)  replace
}


/*----------------------------------------------------------------------
Table A.2: Baseline balance table of parenting knowledge and skills
-----------------------------------------------------------------------*/
{		
orth_out 	belief_able_read_bl belief_able_play_bl belief_read_bl belief_play_bl belief_start_solid ///
			belief_most_imp_food belief_balanced_diet belief_nutr_health belief_nutr_dev belief_nutr_sign belief_nutr_suppl ///
			using  "$outputdir/taba2_balance", by(treatment) covariates(provid*) vce(cluster villid) ///
			pcompare count se numlabel armlabel("Control" "Treatment") title(Table A.2: Baseline parenting knowledge and skills)  replace
}


/*----------------------------------------------------------------------
Table A.3: Baseline balance table of parenting practices
-----------------------------------------------------------------------*/
{			
orth_out 	inv_toy_bl inv_song_bl inv_story_bl inv_book_bl inv_books_bl ch_breastfed_ever initiation continued ch_formula ch_breastfed_duration ch_formula_duration ch_suppl_food dietary_div_bl min_meal_freq_bl ///
			using  "$outputdir/taba3_balance", by(treatment) covariates(provid*) vce(cluster villid) ///
			pcompare count se numlabel armlabel("Control" "Treatment") title(Table A.3: Baseline parenting practices)  replace
}

/*----------------------------------------------------------------------
Table A.4: Baseline balance table for stayers
-----------------------------------------------------------------------*/
{		
orth_out 	cog_score2015 rc_score2015 ec_score2015 fm_score2015 gm_score2015 soemo_score2015 ill_diarrhea_bl ill_tract_bl ill_times_bl hb_bl ///
			ch_anemic_bl height_bl weight_bl waz_bl underweight_bl haz_bl whz_bl ///
			using  "$outputdir/taba4_balance_stayers" if month2016!=., by(treatment) covariates(provid*) vce(cluster villid) ///
			pcompare count se numlabel armlabel("Control" "Treatment") title(Table A.4: Baseline ECD outcomes of non-attritted intention-to-treat population)  replace
}


/*-----------------------------------------------------
Table A.5: ITT effects on psychosocial stimulation
-------------------------------------------------------*/
{

* Panel A: Beliefs
	
	eststo clear 
	
	eststo: reg knowledge_index		treatment 	knowledge_index_bl    		i.provid, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su knowledge_index if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg belief_read			treatment	belief_read_bl 				i.provid, vce(cluster villid) 	
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su belief_read if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg belief_play			treatment	belief_play_bl 				i.provid, vce(cluster villid)
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su belief_play if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg belief_able_read   	treatment	belief_able_read_bl 		i.provid, vce(cluster villid)
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su belief_able_read if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg belief_able_play   	treatment 	belief_able_play_bl    	 	i.provid, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su belief_able_play if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	
					estout using  "$outputdir/tab_ITT_beliefs.xls", replace cells(b(fmt(2))  se(fmt(2))  ci(fmt(2)) ) starlevels(* 0.2 ** 0.1 *** 0.02)  mgroups("Knowledge index" "Knows reading important" "Knows play important" "Able to read" "Able to play", pattern(1 1 1 1 1 1 1 1 1 1 1 1)) collabels(none) numbers("[" "]") mlabels(none) ///
					keep(treatment)	stats(pval1 r2 N cont_mean, labels("P-value" "R-squared" "No. of observations" "Control mean"))											///
					title(ITT Effects on Parental Beliefs) 								///
					note("Expressed in z-scores. Z-scores were obtained by non-parametrically standardizing within age groups. Effects on primary outcomes are estimated intervention effects controlling for baseline scores. Standard errors are adjusted for clustering at the village level and stratification at the township level.")
	
	rwolf 			belief_read belief_play belief_able_read belief_able_play, indepvar(treatment) bl(_bl) strata(provid) controls(i.provid) vce(cluster villid) reps(100) seed(563)

** Panel B: Parental Investment																		
		
	eststo clear 
	
	eststo: reg investment_index  	treatment 	investment_index_bl				i.provid , vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su investment_index if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg inv_toy 			treatment	inv_toy_bl 						i.provid 	, vce(cluster villid) 	
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su inv_toy if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg inv_story   		treatment	inv_story_bl 					i.provid 	, vce(cluster villid)
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su inv_story if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg inv_book   			treatment	inv_book_bl						i.provid 	, vce(cluster villid)
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su inv_book if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg inv_song   			treatment 	inv_song_bl    					i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su inv_song if treatment==0
	estadd scalar cont_mean = `r(mean)'
		
	eststo: reg inv_books   		treatment 	inv_books_bl   					i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su inv_books if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg inv_play_area  		treatment 	inv_play_area_bl   				i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su inv_play_area if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
					estout using  "$outputdir/tab_ITT_investment.xls", replace cells(b(fmt(2))  se(fmt(2))  ci(fmt(2)) ) starlevels(* 0.2 ** 0.1 *** 0.02)  mgroups("Investment index" "Play with toys" "Tell stories" "Read books" "Sing songs" "Play alone" "No. of books" "Play area", pattern(1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)) collabels(none) numbers("[" "]") mlabels(none) ///
					keep(treatment)	stats(pval1 r2 N cont_mean, labels("P-value" "R-squared" "No. of observations" "Control mean"))											///
					title(ITT Effects on Parental Investment) 								///
					note("Expressed in z-scores. Z-scores were obtained by non-parametrically standardizing within age groups. Effects on primary outcomes are estimated intervention effects controlling for baseline scores. Standard errors are adjusted for clustering at the village level and stratification at the township level.")
								
	rwolf 			inv_toy inv_story inv_book inv_song inv_books inv_play_area, indepvar(treatment) bl(_bl) strata(provid) controls( i.provid) vce(cluster villid) seed(563)
	                                             
}	
	
	

/*-----------------------------------------------------
Table A.6: ITT effects on health promotion
-------------------------------------------------------*/
{
* Panel A: Child Dietary Diversity																	


	eststo clear 
	
	eststo: reg dietary_div		 		treatment	dietary_div_bl 				i.provid 	, vce(cluster villid) 	
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su dietary_div if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg iron_rich_index		 	treatment	iron_rich_index_bl 			i.provid 	, vce(cluster villid) 	
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su iron_rich_index if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg food_index		 		treatment	food_index_bl 				i.provid 	, vce(cluster villid) 	
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su food_index if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg ch_staple		   		treatment	ch_staple_bl 				i.provid 	, vce(cluster villid)
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su ch_staple if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg ch_legumes_nuts		   	treatment	ch_legumes_nuts_bl	 		i.provid 	, vce(cluster villid)
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su ch_legumes_nuts if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg ch_dairy		   		treatment 	ch_dairy_bl	    			i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su ch_dairy if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg ch_flesh	  			treatment 	ch_flesh_bl	    			i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su ch_flesh if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg ch_eggs		   			treatment 	ch_eggs_bl	    			i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su ch_eggs if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg ch_vitaminA				treatment 	ch_vitaminA_bl    			i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su ch_vitaminA if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg ch_other_fruit_veg	  treatment 	ch_other_fruit_veg_bl 		i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su ch_other_fruit_veg if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
					estout using  "$outputdir/tab_ITT_food_diversity.xls", replace cells(b(fmt(2))  se(fmt(2))  ci(fmt(2)) ) starlevels(* 0.2 ** 0.1 *** 0.02)  mgroups("Minimum dietary diversity" "Iron rich foods" "Dietary diversity index" "Grains, roots and tubers" "legumes and nuts" "dairy products" "flesh foods" "eggs" "vitamin-A rich fruits or veggies" "other fruits or veggies", pattern(1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)) collabels(none) numbers("[" "]") mlabels(none) ///
					keep(treatment)	stats(pval1 r2 N cont_mean, labels("P-value" "R-squared" "No. of observations" "Control mean"))											///
					title(ITT Effects on Food Diversity) 								///
					note("Expressed in z-scores. Z-scores were obtained by non-parametrically standardizing within age groups. Effects on primary outcomes are estimated intervention effects controlling for baseline scores. Standard errors are adjusted for clustering at the village level and stratification at the township level.")
					
	rwolf			ch_staple ch_legumes_nuts ch_dairy ch_flesh ch_eggs ch_vitaminA ch_other_fruit_veg, indepvar(treatment) strata(provid) bl(_bl) controls(i.provid)  	vce(cluster villid) seed(563) reps(100)																						
	

	
* Panel B: Meal Frequency																			
	
	
	eststo clear 
	
	eststo: reg min_meal_freq	 			treatment 	min_meal_freq_bl    			i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su min_meal_freq if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg meal_freq_index		 		treatment	meal_freq_index_bl				i.provid 	, vce(cluster villid) 	
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su meal_freq_index if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg freq_breastfed		 		treatment	freq_breastfed_bl 				i.provid 	, vce(cluster villid) 	
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su freq_breastfed if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg freq_milk	 	  			treatment 	freq_milk_bl    				i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su freq_milk if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg freq_food		  			treatment 	freq_food_bl    				i.provid 	, vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su freq_food if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
					estout using  "$outputdir/tab_ITT_meal_frequency.xls", replace cells(b(fmt(2))  se(fmt(2))  ci(fmt(2)) ) starlevels(* 0.2 ** 0.1 *** 0.02)  mgroups("Minimum meal frequency" "Meal frequency index" "Still breastfed" "Months breastfed" "Breastfeeding" "Yogurt" "Formula milk" "Milk" "soft, semi-solid or solid food" , pattern(1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)) collabels(none) numbers("[" "]") mlabels(none) ///
					keep(treatment)	stats(pval1 r2 N cont_mean, labels("P-value" "R-squared" "No. of observations" "Control mean"))											///
					title(ITT Effects on Meal Frequency) 								///
					note("Expressed in z-scores. Z-scores were obtained by non-parametrically standardizing within age groups. Effects on primary outcomes are estimated intervention effects controlling for baseline scores. Standard errors are adjusted for clustering at the village level and stratification at the township level.")
					
	rwolf   		freq_breastfed  freq_milk freq_food, indepvar(treatment) strata(provid) controls(i.provid) bl(_bl) vce(cluster villid) seed(563) reps(100)																	
																		
* Panel C: Supplements																		


	eststo clear 	
		
	eststo: reg supplement_index		treatment 	supplement_index_bl 		i.provid  , vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su supplement_index if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg suppl_iron		 		treatment 	suppl_iron_bl	  	 		i.provid  , vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su suppl_iron if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg suppl_zinc		  		treatment 	suppl_zinc_bl	  	 		i.provid  , vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su suppl_zinc if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
	eststo: reg suppl_calc		  		treatment 	suppl_calc_bl	  	 		i.provid  , vce(cluster villid) 
	test _b[treatment]=0
	local sign_treat = sign(_b[treatment])
	display "Ho: coef <= 0  p-value = " ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	estadd scalar pval1 = ttail(r(df_r),`sign_treat'*sqrt(r(F)))
	su suppl_calc if treatment==0
	estadd scalar cont_mean = `r(mean)'
	
					estout using  "$outputdir/tab_ITT_supplements.xls", replace cells(b(fmt(2))  se(fmt(2))  ci(fmt(2)) ) starlevels(* 0.2 ** 0.1 *** 0.02)  mgroups("Supplementation index" "Iron" "Zinc" "Calcium" "Vitamin A/D" "Trace supplements" , pattern(1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 )) collabels(none) numbers("[" "]") mlabels(none) ///
					keep(treatment)	stats(pval1 r2 N cont_mean controls nutrition, labels("P-value" "R-squared" "No. of observations" "Control mean" "Controls" "Nutrition controls"))											///
					title(ITT Effects on Nutritional Supplementation) 								///
					note("Expressed in z-scores. Z-scores were obtained by non-parametrically standardizing within age groups. Effects on primary outcomes are estimated intervention effects controlling for baseline scores. Standard errors are adjusted for clustering at the village level and stratification at the township level.")
					
	rwolf 			suppl_iron suppl_zinc suppl_calc , indepvar(treatment) bl(_bl) strata(provid) controls(i.provid) vce(cluster villid) seed(563) reps(100)																						
}	
	


/*-----------------------------------------------------
Table A.7: Dose-response estimation
-------------------------------------------------------*/
{
              * gen quadratic hh visit var
                gen visit_times_sq = visit_times*visit_times
                
              * generate interactions instruments				
				gen hh_distancekm 	= hh_distance*1000
				gen hh_distancesq 	= hh_distance*hh_distance
				gen dist_treat		= treatment*hh_distancekm
				gen distsq_treat	= hh_distancesq*treatment
				
					lab var visit_times_sq 	"Number of HH Visits^2"
					lab var hh_distancekm 	"Distance from village committee (km)"
					lab var hh_distancesq 	"Distance from village committee (m)^2"
					lab var dist_treat 		"Treatment * distance (km)"
					lab var distsq_treat 	"Treatment * distance (m)^2"
				
		
	* control function estimation 4(distance from village committee)
                reg visit_times treatment hh_distancekm hh_distancesq dist_treat distsq_treat i.provid , vce(cluster villid) 
				test treatment hh_distancekm hh_distancesq dist_treat distsq_treat  
				outreg2 using "$outputdir\first_stage.xls", replace ctitle("HH visits") title("First stage regression results") nocons keep(treatment hh_distancekm hh_distancesq dist_treat distsq_treat ) addstat(F-test, `r(F)', p-val, `r(p)') alpha(0.01, 0.05, 0.10) label addtext(Province FE, YES, Tester FE, YES)
	
                predict v_hat, residuals 
				
				reg visit_times_sq treatment hh_distancekm hh_distancesq dist_treat distsq_treat i.provid , vce(cluster villid) 
				test treatment hh_distancekm hh_distancesq dist_treat distsq_treat  
				outreg2 using "$outputdir\first_stage.xls", append ctitle("HH visits^2") title("First stage regression results") nocons keep(treatment hh_distancekm hh_distancesq dist_treat distsq_treat ) addstat(F-test, `r(F)',p-val, `r(p)') alpha(0.01, 0.05, 0.10) label addtext(Province FE, YES, Tester FE, YES)
	
                predict v_hat_sq, residuals
				 
   
		* Skill Development
	
                eststo  clear
                eststo: reg          skill_dev_index     visit_times                            v_hat              skill_dev_index_bl       i.provid     i.tester_id_2016, vce(cluster villid)  
                outreg2 using  "$outputdir/tab_dose_response.xls", replace  alpha(0.01,0.05,0.1)  label  nocons                 ///                                                                                                        ///
                keep( visit_times   skill_dev_index_bl  v_hat ) title("Dose-Response Relationship")   ctitle("Skill development index score")    addtext(Province FE, YES, Tester FE, YES)
				
				eststo  clear
                eststo: reg          skill_dev_index     visit_times   visit_times_sq           v_hat  v_hat_sq    skill_dev_index_bl        i.provid     i.tester_id_2016, vce(cluster villid) 
                outreg2 using  "$outputdir/tab_dose_response.xls", append  alpha(0.01,0.05,0.1)  label  nocons                 ///                                                                                                        ///
                keep( visit_times   visit_times_sq skill_dev_index_bl   v_hat v_hat_sq) title("Dose-Response Relationship")   ctitle("Skill development index score")   addtext(Province FE, YES, Tester FE, YES)
	
	
		* Health
   
                eststo  clear
                eststo: reg          health_index     visit_times                               v_hat              	health_index_bl        	i.provid , vce(cluster villid)  
                outreg2 using  "$outputdir/tab_dose_response.xls", append  alpha(0.01,0.05,0.1)  label  nocons                 ///                                                                                                        ///
                keep( visit_times   health_index_bl  v_hat) title("Dose-Response Relationship")   ctitle("Health index score")    addtext(Province FE, YES)
				
				eststo  clear
                eststo: reg          health_index     visit_times   visit_times_sq              v_hat  v_hat_sq    	health_index_bl        	i.provid  , vce(cluster villid) 
                outreg2 using  "$outputdir/tab_dose_response.xls", append  alpha(0.01,0.05,0.1)  label  nocons                 ///                                                                                                        ///
                keep( visit_times   visit_times_sq  health_index_bl v_hat v_hat_sq) title("Dose-Response Relationship")   ctitle("Health index score")   addtext(Province FE, YES)
				
			}	
